Podsumowanie analizy danych

W mojej analizie na samym początku przy wczytywaniu danych odrzuciłam atrybuty, które nie są przeze mnie używane w dalszej analizie, ani regresji bądź klasyfikacji.

Następnie z mojego zbioru usunęłam rekordy, w których wystąpiła niewiadoma wartość. Zdecydowałam się na usunięcie, ponieważ zbiór danych jest na tyle duży, że mogłam sobie na to pozwolić zamiast zastępowania tych wartości średnimi, co jest pewnego rodzaju przekłamaniem.

Dodatkowo ze zbioru pozbyłam się tych ligandów, które zostały uwzględnione w zbiorze do usunięcia.

W kolejnym punkcie wykonałam krótkie podsumowanie najważniejszych danych wykorzystywanych w dalszym przetwarzaniu.

Następnie znalazłam 50 najczęściej występujących ligandów i ograniczyłam zbiór właśnie do nich oraz przedstawiłam je na histogramie, posortowane według częstości wystąpienia.

Jednym z wyzwań było przedstawienie korelacji między zmiennymi. Chciałam w tym celu użyć biblioteki corrplot, która przedstawia graficznie macierz korelacji. Niestety przy tak dużej ilości atrybutów musiałam z tego zrezygnować i postanowiłam przedstawić korelacje w postaci tabeli dodatkowo ograniczając wyniki tylko do tych, gdzie korelacja jest większa bądź równa 0.999.

Aby zaprezentować rozkład wartości liczby atomów i elektronów zastosowałam wykresy gęstościowe.

Dziesięć największych niezgodności liczby atomów i elektronów przedstawiłam w formie tabeli, podając dla każdej grupy ligandów uśrednioną różnicę.

Kolejnym z wyzwań było zaprezentowanie rozkładu wartości kolumn part_01. Wynikało to z ilości tych kolumn. W celu drobnej optymalizacji czytelności wykresów zastosowałam krok usunięcia outlier’ów oraz odcięcia prefixu “part_01_” z nazw kolumn, ponieważ występował w każdej z nich. Początkowo planowałam zastosować bibliotekę ggplot2, niestety nie dawała rady nawet przy 10 tysiącach rekordów. Następnie całość przepisałam na bibliotekę plotly, co zadziałało dla 10 tysięcy rekordów, ale niestety dla całości danych wygenerowało raport wielkości ponad 0.4 GB, którego nie byłam w stanie otworzyć w przeglądarce. Ostatecznie byłam zmuszona użyć najprostszego rozwiązania ze statycznymi wykresami.

Regresję liniową wykonałam dla standardowych parametrów. Stratyfikowany podział zbioru na uczący i testowy jako odpowiednio 70% i 30% danych oraz wykorzystanie funkcji lm.

Tworząc klasyfikator zmieniłam stosunek podziału na zbiór uczący i testowy, ponieważ danych było wystarczająco dużo. Dlatego ostatecznie stratyfikowany podział wynosi 99% danych dla zbioru uczącego i 1% dla zbioru testowego. Do uczenia użyłam algorytmu Random Forest z parametrem ntrees = 1000.

Wykorzystane biblioteki

library(dplyr)
library(ggplot2)
library(plotly)
library(tidyr)
library(knitr)
library(caret)
library(data.table)
library(tibble)

Wczytywanie danych

removable_columns <- c("title", "pdb_code", "res_id", "chain_id", "local_res_atom_count", "local_res_atom_non_h_occupancy_sum", "local_res_atom_non_h_electron_occupancy_sum", "local_res_atom_C_count", "local_res_atom_N_count", "local_res_atom_O_count", "local_res_atom_S_count", "dict_atom_C_count", "dict_atom_N_count", "dict_atom_O_count", "dict_atom_S_count", "skeleton_data", "skeleton_cycle_4", "skeleton_diameter", "skeleton_cycle_6", "skeleton_cycle_7", "skeleton_closeness_006_008", "skeleton_closeness_002_004", "skeleton_cycle_3", "skeleton_avg_degree", "skeleton_closeness_004_006", "skeleton_closeness_010_012", "skeleton_closeness_012_014", "skeleton_edges", "skeleton_radius", "skeleton_cycle_8_plus", "skeleton_closeness_020_030", "skeleton_deg_5_plus", "skeleton_closeness_016_018", "skeleton_closeness_008_010", "skeleton_closeness_018_020", "skeleton_average_clustering", "skeleton_closeness_040_050", "skeleton_closeness_014_016", "skeleton_center", "skeleton_closeness_000_002", "skeleton_density", "skeleton_closeness_030_040", "skeleton_deg_4", "skeleton_deg_0", "skeleton_deg_1", "skeleton_deg_2", "skeleton_deg_3", "skeleton_graph_clique_number", "skeleton_nodes", "skeleton_cycles", "skeleton_cycle_5", "skeleton_closeness_050_plus", "skeleton_periphery", "fo_col", "fc_col", "weight_col", "grid_space", "solvent_radius", "solvent_opening_radius", "part_step_FoFc_std_min", "part_step_FoFc_std_max", "part_step_FoFc_std_step")

data <- fread("./all_summary.csv", header = TRUE, drop = removable_columns)

Przetwarzanie brakujących danych

data <- data %>% 
  drop_na()

Usuwanie niepotrzebnych ligandów

deletable_res_name <- c("UNK", "UNX", "UNL", "DUM", "N", "BLOB", "ALA", "ARG", "ASN", "ASP", "CYS", "GLN", "GLU", "GLY", "HIS", "ILE", "LEU", "LYS", "MET", "MSE", "PHE", "PRO", "SEC", "SER", "THR", "TRP", "TYR", "VAL", "DA", "DG", "DT", "DC", "DU", "A", "G", "T", "C", "U", "HOH", "H20", "WAT")
data <- data %>% filter(!res_name %in% deletable_res_name)

Podsumowanie danych

cols_summary <- c("res_name", "local_res_atom_non_h_count", "local_res_atom_non_h_electron_sum", "dict_atom_non_h_count", "dict_atom_non_h_electron_sum")
statistics <- data %>%
  select(one_of(cols_summary))
kable(summary(statistics))
res_name local_res_atom_non_h_count local_res_atom_non_h_electron_sum dict_atom_non_h_count dict_atom_non_h_electron_sum
Length:525666 Min. : 1.00 Min. : 3.0 Min. : 1.00 Min. : 3.0
Class :character 1st Qu.: 4.00 1st Qu.: 30.0 1st Qu.: 4.00 1st Qu.: 30.0
Mode :character Median : 6.00 Median : 48.0 Median : 6.00 Median : 48.0
NA Mean : 13.11 Mean : 98.3 Mean : 13.61 Mean : 101.8
NA 3rd Qu.: 18.00 3rd Qu.: 125.0 3rd Qu.: 19.00 3rd Qu.: 132.0
NA Max. :111.00 Max. :1223.0 Max. :128.00 Max. :1223.0
dim(data)
## [1] 525666    350

50 najpopularniejszych ligandów

popular_ligands <- data %>%
  select(res_name) %>%
  count(res_name, sort = TRUE) %>%
  slice(1:50)

popular_names_vector <- popular_ligands %>%
  pull(res_name)

data <- data %>% filter(res_name %in% popular_names_vector)

Korelacja między zmiennymi

data %>%
  select_if(is.numeric) %>%
  cor() %>%
  as.data.frame() %>%
  rownames_to_column() %>%
  gather(rowname2, value, -rowname) %>%
  filter(value >= 0.9999, rowname != rowname2) %>%
  kable()
rowname rowname2 value
part_00_max local_max 1.0000000
part_01_max local_max 1.0000000
part_02_max local_max 0.9999999
part_00_max_over_std local_max_over_std 1.0000000
part_01_max_over_std local_max_over_std 1.0000000
part_02_max_over_std local_max_over_std 0.9999999
part_00_density_segments_count part_00_shape_segments_count 1.0000000
part_00_shape_segments_count part_00_density_segments_count 1.0000000
part_00_shape_M000 part_00_volume 1.0000000
part_00_density_M000 part_00_electrons 1.0000000
local_max part_00_max 1.0000000
part_01_max part_00_max 1.0000000
part_02_max part_00_max 1.0000000
local_max_over_std part_00_max_over_std 1.0000000
part_01_max_over_std part_00_max_over_std 1.0000000
part_02_max_over_std part_00_max_over_std 1.0000000
part_00_volume part_00_shape_M000 1.0000000
part_00_electrons part_00_density_M000 1.0000000
part_01_density_segments_count part_01_shape_segments_count 1.0000000
part_01_shape_segments_count part_01_density_segments_count 1.0000000
part_01_shape_M000 part_01_volume 1.0000000
part_01_density_M000 part_01_electrons 1.0000000
local_max part_01_max 1.0000000
part_00_max part_01_max 1.0000000
part_02_max part_01_max 1.0000000
local_max_over_std part_01_max_over_std 1.0000000
part_00_max_over_std part_01_max_over_std 1.0000000
part_02_max_over_std part_01_max_over_std 1.0000000
part_01_volume part_01_shape_M000 1.0000000
part_01_electrons part_01_density_M000 1.0000000
part_02_density_segments_count part_02_shape_segments_count 1.0000000
part_02_shape_segments_count part_02_density_segments_count 1.0000000
part_02_shape_M000 part_02_volume 1.0000000
part_02_density_M000 part_02_electrons 1.0000000
local_max part_02_max 0.9999999
part_00_max part_02_max 1.0000000
part_01_max part_02_max 1.0000000
local_max_over_std part_02_max_over_std 0.9999999
part_00_max_over_std part_02_max_over_std 1.0000000
part_01_max_over_std part_02_max_over_std 1.0000000
part_02_volume part_02_shape_M000 1.0000000
part_02_electrons part_02_density_M000 1.0000000

Liczność najpopularniejszych ligandów według nazwy

plot_ligands <- ggplot(popular_ligands, aes(x = reorder(res_name, -n), y = n, fill = n)) +
  geom_bar(stat = "identity") +
  theme(axis.text.x = element_text(angle = 90)) +
  xlab("ligand")+
  ylab("liczność") +
  labs(title = "Liczność ligandów według nazwy")

ggplotly(plot_ligands)

Rozkłady gęstościowe liczb

Atomów

plot_atom <- ggplot(data, aes(x = local_res_atom_non_h_count)) +
  geom_density(alpha = .3, fill = "#00CECB", color = NA) +
  xlab("liczność atomów") +
  ylab("gęstość") +
  labs(title = "Rozkład gęstościowy atomów")

ggplotly(plot_atom)

Elektronów

plot_electron <- ggplot(data, aes(x = local_res_atom_non_h_electron_sum)) +
  geom_density(alpha = .3, fill = "#FF5E5B", color = NA) +
  xlab("liczność elektronów") +
  ylab("gęstość") +
  labs(title = "Rozkład gęstościowy elektronów")

ggplotly(plot_electron)

Największe niezgodności liczby

Atomów

data %>%
  select(res_name, local_res_atom_non_h_count, dict_atom_non_h_count) %>%
  group_by(res_name) %>%
  summarise(atom_inconsistency = mean(abs(local_res_atom_non_h_count - dict_atom_non_h_count))) %>%
  arrange(-atom_inconsistency) %>%
  slice(1:10) %>%
  kable()
res_name atom_inconsistency
CLA 3.5782736
1PE 2.6736000
COA 1.8860182
MLY 1.3122642
NAP 1.2863501
PG4 1.0391591
SEP 1.0007042
NDP 0.9985022
NAG 0.9798955
MAN 0.8681882

Elektronów

data %>%
  select(res_name, local_res_atom_non_h_electron_sum, dict_atom_non_h_electron_sum) %>%
  group_by(res_name) %>%
  summarise(electron_inconsistency = mean(abs(local_res_atom_non_h_electron_sum - dict_atom_non_h_electron_sum))) %>%
  arrange(-electron_inconsistency) %>%
  slice(1:10) %>%
  kable()
res_name electron_inconsistency
CLA 21.634236
1PE 18.233600
COA 14.178825
MLY 9.888050
NAP 8.645994
SEP 8.050000
NAG 7.836841
PG4 7.029679
MAN 6.945505
PLP 6.885058

Rozkład wartości kolumn part_01

remove_outliers <- function(data, na.rm = TRUE, ...) {
  qnt <- quantile(data, probs=c(.25, .75), na.rm = na.rm, ...)
  iqr <- 1.5 * IQR(data, na.rm = na.rm)
  data_no_outliers <- data
  data_no_outliers[data < (qnt[1] - iqr)] <- NA
  data_no_outliers[data > (qnt[2] + iqr)] <- NA
  data_no_outliers[!is.na(data_no_outliers)]
  data_no_outliers
}

plot_part_data <- data %>%
  select(contains("part_01"))

plot_part_data <- plot_part_data %>%
  sapply(remove_outliers) %>%
  as.data.frame() %>%
  drop_na()

names(plot_part_data) <- gsub("part_01_", "", names(plot_part_data))

plot_part_data <- plot_part_data %>%
  gather(part, value, 1:106)

ggplot(plot_part_data, aes(x = part, y = value)) +
  geom_boxplot() +
  facet_wrap(~part, scales = "free", ncol = 6) +
  labs(title = "Rozkład wartości kolumn part_01")

Regresja liniowa

Liczba atomów

data_partition <- data %>%
  select_if(is.numeric)

set.seed(111)
partition <- createDataPartition(
  y = data_partition$local_res_atom_non_h_count,
  p = .7,
  list = FALSE)

data_train <- data_partition %>%
  slice(partition)
data_test <- data_partition %>%
  slice(-partition)

set.seed(111)
fit <- train(local_res_atom_non_h_count ~ ., data = data_train, method = "lm")
fit
## Linear Regression 
## 
## 245042 samples
##    346 predictor
## 
## No pre-processing
## Resampling: Bootstrapped (25 reps) 
## Summary of sample sizes: 245042, 245042, 245042, 245042, 245042, 245042, ... 
## Resampling results:
## 
##   RMSE       Rsquared   MAE       
##   0.1525355  0.9998382  0.03701618
## 
## Tuning parameter 'intercept' was held constant at a value of TRUE
set.seed(111)
prediction <- predict(fit, newdata = data_test)
postResample(pred = prediction, obs = data_test$local_res_atom_non_h_count)
##       RMSE   Rsquared        MAE 
## 0.12337372 0.99991220 0.03643602

Liczba elektronów

data_partition <- data %>%
  select_if(is.numeric)

set.seed(111)
partition <- createDataPartition(
  y = data_partition$local_res_atom_non_h_electron_sum,
  p = .7,
  list = FALSE)

data_train <- data_partition %>%
  slice(partition)
data_test <- data_partition %>%
  slice(-partition)

set.seed(111)
fit <- train(local_res_atom_non_h_electron_sum ~ ., data = data_train, method = "lm")
fit
## Linear Regression 
## 
## 245042 samples
##    346 predictor
## 
## No pre-processing
## Resampling: Bootstrapped (25 reps) 
## Summary of sample sizes: 245042, 245042, 245042, 245042, 245042, 245042, ... 
## Resampling results:
## 
##   RMSE       Rsquared   MAE      
##   0.9670377  0.9998794  0.2496657
## 
## Tuning parameter 'intercept' was held constant at a value of TRUE
set.seed(111)
prediction <- predict(fit, newdata = data_test)
postResample(pred = prediction, obs = data_test$local_res_atom_non_h_electron_sum)
##      RMSE  Rsquared       MAE 
## 0.8638487 0.9999048 0.2498607

Klasyfikator

Przewidywanie wartości res_name

# removable_columns <- c("blob_coverage", "res_coverage", "local_res_atom_non_h_count", "local_res_atom_non_h_electron_sum", "dict_atom_non_h_count", "dict_atom_non_h_electron_sum")
# data_partition <- data %>%
#   select(-removable_columns)
# 
# data_partition$res_name <- as.factor(data_partition$res_name)
# 
# set.seed(111)
# partition <- createDataPartition(
#   y = data_partition$res_name,
#   p = .99,
#   list = FALSE)
# 
# data_train <- data_partition %>%
#   slice(partition)
# data_test <- data_partition %>%
#   slice(-partition)
# 
# set.seed(111)
# fit <- train(
#   res_name ~ .,
#   data = data_train,
#   method = "rf",
#   ntree = 1000,
#   na.action  = na.pass)
# fit
# 
# set.seed(111)
# prediction <- predict(fit, newdata = data_test)
# confusionMatrix(data = prediction, data_test$res_name)